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'^ ' Abstract. 



The nonstationary and steady-state transport through a mesoscopic sample 
CJ , connected to particle reservoirs via time-dependent barriers is investigated within the 

reduced density operator method. The generalized Master equation is solved via the 
^S| ■ Crank-Nicolson algorithm by taking into account the memory kernel which embodies 

^ I the non-Markovian effects that are commonly disregarded. The lead-sample coupling 

takes into account the match between the energy of the incident electrons and the 
f~^ , levels of the isolated sample, as well as their overlap at the contacts. Using a tight- 

^^ ' binding description of the system we investigate the effects induced in the transient 

[■*-. ■ current by the spectral structure of the sample and by the localization properties of 

^D I its eigenfunctions. In strong magnetic fields the transient currents propagate along 

edge states. The behavior of populations and coherences is discussed, as well as their 
connection to the tunneling processes that are relevant for transport. 



PACS numbers: 73.23.Hk, 85.35.Ds, 85.35.Be, 73.21.La 
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1. Introduction and motivation 

The transport properties of semiconductor structures have mostly been studied by 
measuring the steady-state current in response to a constant source-drain voltage 
drop applied on some leads that connect to the sample. This current brings relevant 
information about resonant tunneling processes. The dependence of the transport 
coefficients on various parameters that are varied in the transport experiments (like 
magnetic field, plunger gate voltages, tunneling coefficients) revealed effects that are 
now milestones of nanoscale transport: Coulomb blockade, Aharonov-Bohm oscillations, 
mesoscopic Fano effect etc. 

On the other hand, there is a growing interest on the electron dynamics inside 
quantum dot structures submitted to time-dependent signals applied at the contacts. 
In this case the system exhibits a more complex behavior than in the steady-state regime 
and the quantity of interest is the time-dependent current in the leads. Recent proposals 
for coherent control of electron spin in a quantum dot make use of time-dependent 
signals [1] and the real-time detection of electron tunneling through a quantum dot 
electrostatically coupled to a charge detector has been reported. [2] Another example 
is the pump-and-probe technique proposed by Tarucha et al. in order to extract 
information about the spin relaxation time from transient current measurements. |3] 

From the theoretical point of view the transient current calculations have been 
primarily performed within the non-equilibrium Keldysh formalism. [H O [6l [7] Also, an 
extension of the Lippmann-Schwinger formalism to time-dependent scattering potentials 
was presented in Ref. |8] 

A characteristic feature of the abovementioned experiments is that the system under 
study is in some sense prepared. More precisely, the chemical potentials of the leads are 
such that the first excited state is above the bias window, while the ground state is 
embedded in it. Obviously, the current in the leads cannot capture all the details of the 
dynamics of electrons in the sample. A suitable theoretical description of time-dependent 
transport through mesoscopic systems on which initial conditions are imposed should 
therefore focus on the system itself. 

The natural formal tool to be used is then the density matrix formalism which 
was successfully used in various problems of quantum optics, p] When adapted to 
electronic transport the general strategy goes as follows: i) One starts with several 
disconnected subsystems, i.e. a sample S and some particle reservoirs characterized by 
different chemical potentials; ii) at instant to the sample is coupled to the reservoirs via a 
transfer Hamiltonian Ht which can be in general time-dependent; iii) starting from the 
quantum Liouville equation for the statistical operator W{t) that describes the total 
system one performs a partial trace over the reservoirs and writes down an integro- 
differential equation for the reduced density operator (RDO, defined in Eq. ((HI)). The 
latter is called the generalized Master equation (GME) since it contains both diagonal 
and off-diagonal elements of RDO. We recall that the usual rate equation describes 
the evolution of the diagonal elements of RDO, i.e. the populations. In the GME 
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the effect of the reservoirs on the sample is taken into account through the so called 
memory kernel which contains an infinite sum of time-ordered multiple commutators 
of the type [Ht, [..[Ht,p]] and therefore relates the RDO at instant t to its history at 
previous times. Otherwise stated, in its general form the equation for the RDO is non- 
Markovian. Usually the effect of the leads is taken into account up to the second order 
in Ht which at the physical level describes sequential tunneling processes. The Markov 
approximation assumes correlation functions in the leads rapidly decaying in time. As 
pointed out by Timm[TO] the characteristic time for the decay of correlations in the 
leads is inverse proportional to the applied bias so that in the linear response regime 
the Markov approximation could be again inappropriate. 

The Born-Markov approximation seems reasonable for steady-state calculations of 
the current in the case of a rather large bias but its applicability to transient regime is 
not so clear and has been even questioned recently. [Hj In particular for a rapidly varying 
pulse applied on the leads or at the contacts one cannot assume that the correlation 
functions decay in time. We recall that such a setup is used in experiments with turnstile 
pumps. [H] Also, if one computes higher moments (e.g. noise) non-Markovian effects need 
to be included. [151 [IB] 

In view of these considerations the aim of this paper is on one hand to investigate 
the time-dependent transport in mesoscopic structures by solving the GME without 
using the Markov approximation. On the other hand we propose an implementation 
of the generalized master equation which allows us to take into account the geometry 
of the sample and uses its spectral properties in order to set an "effective" size for the 
reduced density matrix to be computed numerically. We believe that this represents an 
important step forward because it opens the way to study larger systems and capture 
geometrical effects and details about the electron dynamics inside the system itself. 

In order to set the general framework we give below a brief survey of several versions 
of the RDO method that have been proposed in the context of quantum transport. In 
contrast to the quantum optics where the reservoir is a bosonic environment describing 
the radiation field, in transport problems the system is coupled to particle reservoirs. 
Bruder and Schoeller [17] established a quantum Master equation for the diagonal 
elements of the statistical operator by performing a systematic perturbative expansion 
in powers of Ht- Each term in this expansion corresponds to a tunneling process. 
Their calculations were primarily focused on the steady state regime and emphasized the 
interplay of sequential tunneling and inelastic cotunneling processes (the latter ones are 
described by taking into account fourth order terms in the transfer Hamiltonian) . In the 
real-time diagrammatic approach [18] one writes the Hamiltonian of the central region in 
terms of many-body states that have to be computed in the presence of electron-electron 
interaction. For few-level systems one can therefore use exact diagonalization techniques 
[T9] in order to investigate elastic or inelastic cotunneling or, as is done in Ref. [17], to 
describe the Coulomb interaction within the orthodox model. This is somehow different 
from the standard perturbative calculations within the nonequilibrium Green-Keldysh 
formalism [20] where the interaction term is included as a two-particle operator in the 
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Hainiltonian and therefore an interaction self-energy has to be computed. 

Later on Gurvitz and Prager [21] reahzed that in the hmit of a high bias the density 
matrix obeys modified rate equations resembhng the Bloch equations. A common 
feature of the rate equation approach is that the level broadening due to the coupling 
to the contacts should be included 'by hand' in the equations when integrating over 
energy. In a recent paper [22] Tokura et al. investigated interference effects in parallel 
quantum dot systems in steady-state regime using both the Keldysh approach and the 
Bloch equations. Pedersen and Wacker [23] developed a scheme that holds for arbitrary 
bias and goes beyond the rate equation approach. The matrix elements of the statistical 
operator are computed within the Markov approximation in the wide-band limit. The 
authors find out that in the steady-state regime the Born-Markov approximation for 
the RDO and the non-equilibrium Green's function formalism (NEGF) lead to similar 
currents, while in the transient regime the two approaches give different results. 

Another implementation of the RDO method for the electronic transport in 
quantum dot systems was presented by Harbola et al. [21] Using the Born-Markov 
approximation and the wide-band limit the authors have computed both the Fock space 
populations (FSP) and coherences (FSC) for a two-level noninteracting quantum dot. 
The Master equation that is solved in their case describes the projection of the density 
operator on the n-particle sector of the Fock space; this procedure was introduced 
by Rammer et al. [25] in the context of quantum measurement theory. It was also 
shown that within the rotating wave approximation the coherences are decoupled from 
the populations and then the effects of the former on the steady-state currents is not 
included. A thorough comparative analysis of various GME methods was presented 
by C. Timm. [TO] The master equation derived by Schoeller et al. within the real-time 
diagramatic methods [18] is shown to coincide to the Wangsness- Bloch version of the 
GME. 

In a recent work Vaz et al. [H] used Laplace transform methods in order to compute 
the Redfield tensor that characterizes the memory kernel and presented numerical 
simulations for the Fock space coherences in the non-Markovian case. The main 
statement of this work is the existence of long-lived FSC even in the steady-state. 
Finally one should emphasize that the GME method was also employed for studying 
the effect of a periodically oscillating signal applied on the sample itself rather than on 
the leads. [T2l [13] In the present work the time-dependence appears only in the transfer 
Hamiltonian describing the contacts between the leads and the sample. 

The paper is organized as follows. Section II describes the model and presents the 
main equations, some of the formal details being given in Appendix A. In Section III we 
present the applications of the method and discuss the numerical results. Conclusions 
are summarized in Section IV. 
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Figure 1. A schematic picture of the system we have studied. Two (semi-infinite) 
one-dimensional leads are attached to a two-dimensional sample described by a A^ x M 
sites. The solid dots indicate the contact sites; each site i of the sample is characterized 
by two indices (xi, jji) where Xi ~ 1, .., N and jji — 1, .., M. 



2. Theory 

2.1. The model Hamiltonian 

We consider a mesoscopic sample that is coupled to two leads (particle reservoirs) at 
the initial instant to = 0, but decoupled at earlier time. The reservoirs have different 
chemical potentials. We have therefore three subsystems: the two semi-infinite leads 
I = L,R (Left and Right), and the central sample S. The transport problem concerns 
the evolution of this open quantum system for t > to- 

We denote by ipl^ and e\q) the eigenfunctions and eigenvalues of the single particle 
Hamiltonians hi and hji describing the semi-infinite leads. The sample Hamiltonian 
hs has eigenfunctions 0„ and eigenvalues En- The single particle Hamiltonian of the 
disconnected system is ho = hi + hR-\- hs. We use small letters for these Hamiltonians 
in order to distinguish them form their second quantized form which we shall denote 
below by capital letters. 

Our method can be implemented both for continuous or discrete models. In this 
work we shall present only the tight-binding case. The central region is therefore 
described by a two dimensional lattice Hamiltonian and the leads are modeled as one 
dimensional tight-binding chains. The Hamiltonian of the sample in the coordinate 
representation reads: 

/^5 = $^A,K)(2| + 5^(t5e^^-K)(j| + /i.c), (1) 

where i,j are nearest neighbor sites in the sample, Aj are on-site energies and the Peierls 
phase attached to the hopping parameter is describes a constant perpendicular magnetic 
field. 

In order to describe the coupling between the two subsystems we shall add a 
perturbation to ho. We start from the well known single-particle form of the transfer 
Hamiltonian: 

hT=Yl Xi{t)Vi{\Oi){zi\ + h.c), (2) 

l=L,R 
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where 0; is the site of the lead / which couples to the contact site ii in the sample (note 
that the index ii is identified once we established a labelling of the sites in the central 
region). The time- dependent coupling to the leads is characterized by the switching 
functions x'- Note that in general x^ 7^ X^ for t > t^ and we are not restricted to the 
sudden coupling of the leads as is done in previous works. From the physical point of 
view time-dependent couplings can be realized by applied radio- frequency signals to the 
metallic gates between the sample and the leads, as is done for example in the turnstile 
pump experiments. [26] The constants Vi represent the coupling strength to the Z-th 
lead. 

Since we are dealing here with an open system with variable number of particles 
(recall that the semi-infinite leads simulate particle reservoirs) it is mandatory to switch 
to a many-particle Hamiltonian, although in the present work we shall completely neglect 
the Coulomb interaction (we discuss this point further in Section IV). According to the 
general rules of second quantization [27] a basis in the Fock space JF of the coupled 
system can be constructed starting from the eigenfunctions ipq and 0n- One defines 
creation and destruction operators for electrons in the leads c^^ (c^/) and in the sample 
d\^ (dn)- Then the second- quantized total Hamiltonian reads as follows: 

^W= Yl dqe\q)c\fqi + ^Endidn 

+ E E I dqx\t){Tlc\idn + h.c), (3) 

l=L,R n -^ 

where the coefficients T^'^ are given by: 

7^^n = ^<(O)0n(^O, (4) 

The eigenfunctions of the sample are numerically computed while the wave functions 
^/'g are known analytically (note that they are real and do not depend on the lead index 
/ as we take identical leads): 

,, , , sin(o(m + 1)) „ , . 

V' m =^==^=^, e, = 2tLCOsq. 5 

V2tL sm q 

In the above equation t^ is the hopping energy of the leads. Of course one could 
consider more complicated couplings, taking into account more sites from the central 
region coupled to two-dimensional leads. A similar way of constructing coupling matrix 
elements depending on junction configuration was proposed in Ref. [3H]. The integral 
over q counts the momenta of the incident electrons such that e\q) scans the continuous 
spectrum of the semi-infinite leads. 

The third term in Eq. ([3]) is the so called transfer or tunneling Hamiltonian. It 
has been introduced in the early days of electronic quantum transport and thoroughly 
discussed in a series of papers. [28l EHl EU [30] The tunneling Hamiltonian was extensively 
used within the non-equilibrium Green-Keldysh transport formalism. Usually the wide- 
band limit approximation is assumed and then the energy dependence of the coupling 
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coefficients is neglected. [20] More important details of the contacts like position or 
width are also omitted. 

In the present approach the coefficients r^„ computed from Eq. (jlj) contain three 
features: i) The dependence on energies e'' and E^ (through ipl^ and 0^). ii) The precise 
location of the contacts between the leads and the sample (i.e. the sites ii. iii) The 
probabilities \4'g\'^ and |0„p to have electrons at the contact sites. 

2.2. The generalized Master equation (GME) 

Having introduced the second quantized Hamiltonian H{t) we now define the statistical 
operator of the open quantum system as the solution of the Liouville equation: 



iW{t) = [H{t), Wit)], Wit < to) = PLPRPS, (6) 



where 



^-(3{Hi-tiiNi) 

P^ " Tr/{e-'3(^'-'^'^')}' ^ ' 

In the above equation ps is the density operator of the isolated system (that is, for times 
t < to) and serves as an initial condition for the RDO. ni and Ni denote the chemical 
potential and the occupation number operator of the lead /, pi being the equilibrium 
statistical operator of the disconnected lead /. The trace at the denominator is taken 
in the Fock space of the leads. The RDO is defined as the (partial) trace on the Fock 
space of the leads: 

pit) = TTLTTnWit), pito) = Ps. (8) 

The main problem is to find, under suitable approximations, the matrix elements of 
pit) with respect to a basis in the Fock space J-'s of the sample. One way to deal with 
this problem is to compute conditional reduced operators acting in different n-particle 
sectors of the Fock space (see for example [35] ). Moreover, Li et al. [35] proposed a 
factorization for the full density matrix ipit) = YlnP^"'^ ® Pleads) which generalizes the 
usual Born-Markov approximation. [23] In the present approach we do not impose an 
equilibrium state on the leads after the coupling is switched on, which would mean to 
take Wit) = pip^pit). While the steady-state currents are most likely not affected by 
this ansatz the transients are expected to be different when computed within the two 
approaches. 

We shall use the occupation number basis constructed from the single-particle states 
{4>i} of the isolated system. Then a many-body state u reads as: 

|z/) = |^r,^^,..,C-), (9) 

where the number i'^ indicates if the n-th single particle state is occupied (ij^ = 1) or 
empty (zj^ = 0). The corresponding energy of the many body state is denoted by 8,^ and 
is given by the sum of the occupied single-particle levels, i.e Si, = X^n-^^i^n- 

It is clear that the size of the reduced density matrix becomes already very large 
if the central region accommodates A^ ~ 20 electrons and for A^ ~ 50 it seems quite 
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Figure 2. The initial configuration in the many- level quantum dot for a given pair of 
chemical potentials on the leads fJ-LifJ-R. and a 'gap' 6. The occupied levels are marked 
by thick lines. 



impossible to compute the entire matrix, even within the Markov approximation. On the 
other hand one can easily accept that the number of the many-body states (MBS) that 
are relevant to the transport problem is actually much smaller, and at low temperatures 
is controlled by the bias applied on the leads. In the present model the bias is included 
as the difference between the chemical potentials of the leads i.e. eV = fJ^i — f^R, a 
procedure which is also used in the Keldysh formulation of electronic transport. [H [6] 

Suppose now that at instant to the density operator of the central region is such 
that the first A^o single-particle states are occupied and all the higher states are empty, 
that is: 

p(to) = Wo){i^o\, Wo) = |M^,o,o, ) (10) 

A^o states 

where z/q is just the label of the selected many-body state. Moreover, let us consider that 
the bias window is fixed such (ir — En^^ > 6 ioi a positive 5. When the leads are plugged 
to the central region the following scenario is expected: i) The lowest Nq levels remain 
occupied and will not contribute to transport as long as the frequency of the coupling 
signal x' is small compared to the gap 6; ii) It is reasonable to assume that electrons 
tunnel through the dot only via the levels located in the energy range [fiR — 6, fii + S]; 
iii) In the transient regime the occupation numbers of these states will depend on time 
and will eventually settle down in the steady-state regime. Given this setup it is clear 
that there are only (A^max — Nq) single-particle states which are active in the transport 
process and consequently it is sufficient to compute only the matrix elements of the 
RDO for the A/" = 2^™='''^^" many-body states having the following form: 

|z/) = |M^,^J^O + l'--^5^n.ax,0,0, ) (11) 

A^o states 

Let us mention that another interesting initial condition for the density operator is 
p{to) = |i^ex)(t'ex| where: 

|z/ex) = 1 1,1,....1,0,0, ..,im,0,...), im = l (12) 

A^O states 

namely the one in which besides the lowest Nq occupied levels there is an electron on 
a higher level Ej (i.e. the initial states is excited). The decay of the state |z/ex) as the 
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coupling to the leads evolves could be related to the pump-and-probe experiments in 
Ref. [3]. 

We proceed now with the equation of motion for the RDO. It is useful to 
introduce the notation t/o(t, s) = e"**^*"*)^" for the unitary propagator associated to 
the disconnected system {Hq = Hs + Hl + Hr). Using the superoperator method 
developed by Haake [36] we end up with the following GME for the reduced density 
operator up to second order in the tunneling Hamiltonian (we give more details of the 
derivation in Appendix A): 

m = -{[Hs,p{t)] 

-^TiLTiRilHTit), j ds[Uo{t, s)HT{s)Uo{t, s)\ Uo{t, s)p{s)Uo{t, s) Wk]]} (13) 

Note that Uq acts on the entire Fock space and cannot therefore be permuted as a whole 
inside the partial trace; nevertheless, one can do so for e~''^^~^'^^^~^^^' . As a next step 
we rewrite Ht and Hs in terms of many-body states and then work out the double 
commutator in Eq. flT3|) . Using the completeness relation Yla l'^)('^l = 1 we have: 

HAt) = E E [dqx'mT^M^)(f^\^^ + h-C.), (14) 

l=L,R a,/3 -^ 

where we have introduced a scattering operator T acting in the Fock space of the system: 

'^U'i)=J2'^n,WU^)- (16) 

n 

It is clear that T^jji^q) describes the 'absorption' of electrons from the leads to the system 
and changes the many-body states of the latter from j3 ^>- a. Note that in the numerical 
implementation the index n counts only those single-particle states within the active 
energy interval and that in order to have a nonvanishing T the number of electrons in 
the many-body states a and /? have to differ by one. 

Replacing (fTSll in (fT3|) and using the well known identities: 

e^*^'c,,e-^*^'=c,,e-^^'('')*:=c,,(t), 

e^tH,^t^^-..H, ^ ^t^^.H.)* ;= gt^(t), (17) 

as well as the correlation functions of the leads {fi{B\q)) denotes the Fermi function 
that characterizes the lead /): 

T^rdpMtHiit')} = e-^(*-*')^'('^)5(g - k){l - Me\q))) 

TT,{prcl,{t)dM{t')} = e^(*-*')^'(^)5(g - k)Me\q)) 

one writes the GME into a compact form: 

m= -'^[Hs,pit)] 

-T2T. fdqx'm[Ti,QAt)] + h.c) (18) 
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where we have introduced two operators: 

J to 

U,iis) = e-^^ (7;V(^)(1 - /') - Pis)r,^f) e-""' 

For the siniphcity of writing in the above equation we omit to write the energy 
dependence of the Fermi function while keeping only the lead index /. Equation fllSp is 
the main formal result of the paper. It leads to a finite system of A/" coupled integro- 
differential equations for the matrix elements {a\p{t)\j3) of the RDO. Note that the 
commutator structure leads to the conservation of the trace over the sample states, i.e. 
TTsp{t) = 0. All the tunneling processes of second order in the transfer Hamiltonian are 
included in Eq. ( TT8l) and one can identify loss and gain terms contributing to a given 
matrix element of the RDO. Remark that both elastic and inelastic tunneling are taken 
into account. Electrons in a given state of the sample are allowed to tunnel out in a 
state q of the lead at time t and to tunnel back in a different state at time t. Another 
important feature is that the sign change of the coupling matrix elements for subsequent 
levels is fully taken into account in the above GME. For example one can easily check 
in the Eq. flT3|) that terms like T/pT^'/' or T^ pTiil — /') contain products of the form 
T^mqi'^nq)* ^^at cau havc different signs. Such terms describe processes in which one 
electron enters the dot on the m-th level while another one leaves the dot from the n-th 
level. Their role in transport was recently emphasized by Amir et al. [37] 

At this point one can take further approximations on the GME in order to put 
it into a Lindblad form (see for example Ref. ([23])). First one applies the Markov 
approximation in which p(s) is approximated by e*(*~*)'^op(t)e~**^*~*-''^°, the time integral 
is extended to infinity and calculated via a principal value formula. Then one can take 
either the rotating wave approximation or the limit of high bias window (i.e. fL = l and 
Ir = 0). From the physical point of view this limiting case means that the electrons 
can neither flow from the sample to the left lead nor from the right lead to the sample. 
In the present work none of these approximations are needed. 

Once we have the RDO it is possible to compute the statistical average of the charge 
operator in the coupled sample Qs = e ^„ dl^dn- 

{Qsit)) = Tr{Wit)Qs} = TTs{pit)Qs} 



the traces being now assumed in the Fock space. Similarly one introduces the charges 

Ml 

dt 



Ql,r in the leads. We define the net currents in the leads as follows: Jl(^) = — ^ and 



Jnii) = ~M- ^'^ therefore have J/, > if the electrons flow from the left lead towards 
the sample and Jr > if they flow from the sample towards the right lead. In the 
transient regime the sign of the net currents can change. The continuity equation reads 
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= EE<<HpWI^)- (20) 

n u 

Using the GME, Eq. (TT8|) . one can easily identify the contribution of each level n to the 
currents in the left and right lead: 



■Ji,n = -^J2''nJdqx\t)HTi,n,iit)] + h.c.\u) 



(21) 

where zj^ = 0, 1 specifies if the n-th single particle state is occupied or empty. Observe 
that the currents are expressed only in terms of the diagonal elements of p but this does 
not exclude contributions from the off-diagonal elements as well, because all matrix 
elements are coupled in the GME. In order to solve the GME numerically we use the 
Crank-Nicholson method. [3Bj The time is discretized and the first derivative of the 
density operator with respect to time is evaluated as the mean value of the forward (or 
right) and backward (or left) derivatives at the same time point. The time step was 
chosen (and tested) to be sufficiently small in order to capture all physical details of 
the time evolution of the density operator. On the RHS of Eq. flTSj) we approximate 
p{tk+i) by p(tfc) and then perform iterations until a convergence test is fulfilled. The 
time integration included in the operator Vtqi is done recursively: once we know p(tj) for 
any i < k + lwe update the integral by adding the value of the integrand corresponding 
to pitk+i- At any step of the iteration we check the conservation of probability. It 
is known that the form of the GME that we use here, in the lowest coupling order, 
does not guarantee the positivity of the diagonal elements of p, i.e. the probabilities of 
many-body states. In the numerical simulations presented in Section III we checked the 
positivity at each time step. We find that by increasing the coupling constants Vi some 
of the populations could take slightly negative values especially in the transient regime. 
In contrast, for a given sample the variation of the bias around the active region or a 
very fast switching of the coupling do not damage the positivity. 

3. Numerical simulations 

The first sample model is a 5 x 10 lattice which is large enough to exhibit the well 
known Hofstadter spectrum (see Fig 3 (a)) when a strong perpendicular magnetic field 
is applied. [40j The Dirichlet boundary conditions for the two-dimensional discrete 
Laplacian lead to the formation of edge states. [H] The leads are attached at diagonally 
opposite corners of the sample. The magnetic fiux is $ = 0.2 and we take pl = —1.5, 
Pr = —2.25. We consider four active states for transport denoted by Ei,..E4^ and 
distributed as follows: two states are located within the bias window, one state above 
and one state below the window. The rest of the spectrum is separated from this group 
of states by a gap AE ~ 0.5. Therefore the transport properties will be computed from 
a reduced density matrix which accounts for 16 many-body states. The temperature 
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Figure 3. (Color online) (a) A part of the Hofstadter spectrum; the two horizontal 
lines mark the values of the chemical potentials in the leads, i.e. fiL — —1-5 and 
/ifl = —2.25. (b) The effective coupling between a state n in the relevant window and 
a state in the leads having energy Eq. The maxima correspond to resonant tunneling, 
that is En ~ Sq. Other parameters: Vl = Vr, = 1, $ = 0.2. 



is very small, kT = 10~^. The switching functions Xl,r are identical, and describe a 
smooth coupling to the leads XLfiif) = (1 — ^-if^^ )- Note that the parameter 7 decides 
how fast we establish the coupling between the two subsystems (7 = 1 if not stated 
otherwise). As for the coupling strength we used Vl = Vr = 1. The important coupling 
parameters are actually T^^ and one can see in Fig. 3(b) that each state is differently 
coupled to the contact sites. From Eq. (I5l). The shape of |Tg„p as a function of Eg is 
easily see to be like A/4t| — E'^ 

In order to discuss the properties of the RDO we introduce a labeling of many-body 
states. It is clear that the only occupation numbers that will be changed when the sample 
is coupled to the leads are those associated to the active levels. The occupation numbers 
of the active single-particle states define the many-body state |z/), and we shall omit the 
occupation numbers of the frozen (non-active) states. With four active states we denote 
|1) = |0000),|2) = |1000),|3) = |0100),|4) = |1100),|5) = |0010), |6) = |1010),|7) = 
|1100), |8) = |1110), |9) = |0001), etc. (consecutive binary numbers written from right 
to left). 

The interpretation of the sign of the current is: if JL,n and Jr,„ are positive the 
charge flows from the left lead towards the sample and from the sample towards the 
right lead. In the steady state we obtain JL^n — JR,n = for any n. The currents are 
given in ets/h units and the time is expressed in h/ts where ts is the hopping energy in 
the central region; we also take ts as the energy unit. If one considers an effective lattice 
constant a = 10 nm and the effective mass of GaAs in the definition of the hopping 
constant ts = fv^ /2m*a^, it turns out that the time unit is 1 ps and the current unit is 
20 nA. The hopping energy on leads ti = 2tsi in order to match the spectrum of the 
ID lead (i.e [— 2ti, 2ti\) to the spectrum of the central region. 

The active region that we consider first contains only the middle 4 states in Fig. 
3(b). The remaining two were included in order to to check later on that by taking 
two more states in the active region (one below the bias window and one above it) the 
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numerical results are not altered. The coupling depends on the amplitude of the wave 
function. It is clear for example that a state may not contribute to the current even if 
it is energetically inside the bias window if the wave function vanishes at the contact. 
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Figure 4. (Color online) (a) The currents transmitted through the states below the 
bias window (2nd and 3rd), (b) The currents associated to the two states below and 
above the bias window (1st and 4th). (c) The average charge on each level within the 
active region. Other parameters V = 0.75, Vl — Vr — 1. 



First we assume that before the coupling to the leads was established the four 
levels in the active region were empty, that is p(0) = |1)(1|. The details of the electron 
dynamics can be extracted from the currents associated to each level. In Fig. 4(a) we 
show and compare the transients in both leads associated to the two levels within the 
bias window (see the curves corresponding to the labels in upper right corner of the 
figure). In view of further analysis we include as well the currents in the right lead 
when the left lead is disconnected (the two curves corresponding to the labels in the 
lower right corner of the figure). Electrons from both leads can tunnel from or into 
these states and the difference between the chemical potentials leads eventually to equal 
currents in the steady state. In the transient regime the currents in the two leads 
behave however differently: Jl increases abruptly at short times with a bigger slope 
for the level whose coupling to the contact is stronger, while the current flowing into 
the right lead is delayed. Moreover, this delay depends on the state which carries the 
current. The 3rd level starts to transmit charge earlier (for t ~ 5) while the 2nd needs 
more time to inject electrons into the right lead (for t ~ 20). Fig. 4(b) gives the currents 
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passing through the 1st and 4th level, which are outside the bias window. The following 
features are noticed: i) The lowest level absorbs charge from both leads (hence J^i > 
and Jr^i < 0). ii) The currents decrease slowly to zero giving no contribution to the 
steady-state current; iii) The current of the 4th level which is located slightly above the 
bias window oscillates with both positive and negative values because in the transient 
regime this level can gain or loose charge as well; note however that the amplitude of 
the oscillation decreases in time. 

One can also see that the steady state regime is reached faster by the two states 
located within the bias window. This can be better seen in the occupation number of 
the corresponding levels which we present in Fig. 4(c). There are basically two regimes 
for the two levels in the bias window. 
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Figure 5. (Color online) The current transmitted in the right lead through various 
states in the case of an active region containing six levels. Other parameters Vl = 
Vu = 1, kT — 10^"'. The current carried by the 6-th state is vanishingly small and is 
not given in the figure. 



As the system opens they are charging at first by absorbing electrons mainly from 
the left lead (note that the right lead provides as well a small amount of charge for a 
short time); later on the net current in the right lead becomes positive and the steady 
state corresponds to a constant occupation number for each level. This behavior of the 
occupation numbers is consistent with the numerical calculations we performed recently 
using the Green-Keldysh formalism. [6] 

The delay of the current in the right lead which we noticed above could only be 
associated to the time needed for the electrons to propagate along the system and tunnel 
through the contacts. Indeed, when the coupling to the left lead is switched off the time 
needed to get a positive (i.e. outgoing) current in the right lead is almost the same as in 
the two-lead geometry. When Vl = electrons enter into the sample and spend some 
time there before being expelled in the same lead. Nevertheless, in the steady state no 
current is generated in the single lead geometry, and the transient oscillations vanish. 

As for the lowest level one can see that in the steady state it is almost full but the 
time needed to achieve the maximum filling exceed by far the time needed to have the 
2nd and 3rd levels half-filled. This behavior is expected because the lowest level has the 



Geometrical effects and signal delay in time- dependent transport at the nanoscale 15 

smallest coupling to the leads so that the tunneling processes are slower. We shall see 
below that plugging the lead to a contact site where jT^ip is larger will accelerate the 
filling process. Obviously the total current will not reach the steady-state unless each 
partial current does, even if the corresponding active state is outside the bias window. 
Now we take two more single particle states in the active region and compute the 
transients by taking the initial state p(0) = |2)(2| where |2) = |100000) (see Fig. 5). 
This choice allows a comparison with results given in Fig. 4. The 1st and the 6th states 
carry a vanishingly small current while the remaining currents are quite similar to the 
ones in Fig. 4, which justifies the restriction to 4 single particle states in the active region. 
This results could be expected because the couplings Tgi and T^g are the smallest one 
(see Fig 3(b)). 
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Figure 6. (Color online) (a) The total current current in the right lead as a function of 
time for different initial configurations of the isolated sample: solid line - /3(0) — |1)(1|, 
long-dashed Hne - /9(0) = |4)(4|, dotted hue /o(0) = |9)(9|. (b) The currents associated 
to the two states below and above the bias window when the initial state of isolated 
system is given by jo(0) = |9)(9|. Other parameters Vz, = Vr = 1, kT — 10^''. 



The next step of our study is to look at the transients obtained when the levels from 
the active interval are already occupied. Here we shall use the main advantage of the 
reduced density matrix method, that is, to take different initial states for the sample. 
We remind here that in the non-equilibrium Greens' function formalism the only possible 
initial state of the disconnected sample is the vacuum; other initial configurations are 
not naturally implemented. In Fig. 6(a) we compare the total currents in the right lead 
obtained for three initial configurations: p(0) = |1)(1| (empty system), p(0) = |4)(4| 
(the lowest two levels completely occupied) and p(0) = |9)(9| (the highest level from 
the active region occupied). One notices at once that the steady-state currents do not 
depend on the initial configuration but the transients behave differently. First of all, 
when we start with occupied states in the sample there is no delay in the onset of a 
positive current in the right lead, because the electrons immediately tunnel to the leads; 
actually, when p(0) = |9)(9| there is a current flowing to the left lead as well because 
the fourth level is located above the bias window. Secondly, the steady state is achieved 
faster for the conflguration |4) because the lowest level is fllled and remains so up to very 
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small oscillations (not shown). We have checked that the currents Jr^3 and J^4 coincide 
with the ones associated to the initial configuration |1). The relaxation time of the state 
|0001) can be traced back from Fig. 6(b) which shows that the currents Jr^4 and Jl,4 
vanish at t ~ 100; the corresponding occupation number iV4 vanishes also there. Since 
the relaxation process coexists with the filling process of the level located within the 
bias window a local minima appears in the total current followed by a smooth increase 
towards the steady-state value. The turning point in J^ means that the currents fiowing 
to the right lead through the bias window compensates the decrease in the relaxation 
current Jr^4. More importantly, the decay of this state increases the occupation of the 
lowest level. This can be seen in Fig. 6(b) where the currents entering the lowest level 
for the configuration |9) is also. Comparing with the same currents in Fig. 4 it is obvious 
that in this case there is more charge entering the 1st level and that the currents vanish 
faster. 

We discuss now the behavior of the reduced density matrix. The diagonal elements 
refiect the probabilities for having certain occupation numbers in the single particle 
active states. It is clear that some of these configurations are more likely to be realized. 
Figs. 7(a) and (b) show the evolution of some diagonal matrix elements of the RDO. 
Again, we start with the sample in the vacuum state, p{to) = pu, where pu is the 
probability to have all the levels from the active region empty. It decreases from 
the initial value to zero because the levels are populated as the coupling to the leads 
strengthen; the numerical data suggests an exponential decay. 

One can identify several regimes for the populations. At short times the most 
favorable nontrivial populations contain just one electron. Otherwise stated, the system 
spans the single-particle sector of the Fock space. Physically this means that two 
levels cannot be simultaneously occupied shortly after the coupling is switched on. The 
occupation probability of levels with the lowest energies is expected to be higher, i.e. 
P33 > P55- The lowest levels absorb charge from both leads, while the higher ones are 
populated mostly due to tunneling from the left lead. This is true in our case only 
after a short time, but because of the energy dependent coupling we are using, after a 
longer time we see the opposite result. In a second regime the system is more likely to 
be in a state from the two-particle sector of the Fock space, namely |4) = |1100) and 
|6) = 1 1010) (Fig. 7(b)). The corresponding diagonal elements P44 and pee increase with 
time. At even later times the system can accommodate three electrons on the first levels 
so P88 increases as well. In the steady-state regime four configurations have a significant 
probability: |1100), |1010), |1000) and |1110). It is obvious that by pairing these states 
one recovers the tunneling processes that contribute to the steady-state currents in the 
leads. For example, switching from |1100) -^ |1000) implies that one electron tunnels 
out from the 2nd level in the active region. 

The non-Markovian nature of the system implies that off-diagonal elements (the so 
called coherences) of the RDO could develop in time, even if we start from a diagonal 
density operator at t = 0. A nonvanishing coherence Pa/sit) means that at instant t 
the state of the system cannot be completely described by occupation probabihties. 
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Figure 7. (Color online) Time-evolution of the most relevant populations for the 5x10 
sites system, (a) Single-particle configurations decay in the long time limit, (b) The 
populations that do not vanish in the steady state (see the text). The other parameters 
are as in Fig. 6. 



The coherences can be shown to vanish in the steady-state if, on top of the Born- 
Markov approximation, one also takes the rotating wave approximation. [24] This is not 
the case in the GME case and moreover, the diagonal and off-diagonal elements of 
the density operator are coupled. Note however that the second order approximation 
with respect to the transfer Hamiltonian that we keep here implies that coherences 
between many-body states with different particle numbers are excluded. This is because 
TrjpiCgjcJ;} = TT{pLCqiCqi} = 0. Wc show in Figs. 8(a) and (b) the behavior of the 
nonvanishing coherences which are complex quantities and satisfy the relation p*^ = p^a- 
In general both the imaginary and real parts have oscillations. Nevertheless, some 
matrix elements p^a vanish in the steady state (see Fig. 8(b)) while some settle down to 
a non- vanishing value and contribute indirectly to the current, via the diagonal elements 
Paa- We notice also that coherences do not appear simultaneously and their oscillations 
behave differently. For instance, p23 starts earlier than p^Q and its oscillations decrease 
in time; in contrast, although p^Q exhibits mild oscillations at short times it gradually 
increases and clearly exceeds p23 in the steady state. The explanation of this behavior 
lies in the dynamics of the occupation numbers in the sample. On one hand, the state 
|3) = |0100) is realized with probability 15% at short times only (see Fig. 7(a)) but 
then this state become less probable. This is why p23 decreases at t ~ 125. On the 
other hand both states that appear in P35 have higher probability (~ 20%) in the long 
time limit. Similar arguments explain why the other two coherences shown in Fig. 8(b) 
vanish in the steady state: they imply sequences of occupation numbers that are not 
expected in the long time limit. Moreover, since the probability to find the system in 
the state |9) decreases much faster that the one associated to the state |3) it is clear 
that P95 decays well before P35. 

The definition of the coupling coefficients T^^^ (see Eq. (jlj)) implies that by choosing 
different contact regions one could obtain different currents as for any q, n the associated 
eigenf unctions may overlap differently. In Fig. 9(a) we show |Tqip and |Tg3p for 
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Figure 8. (Color online) Time-evolution of the most relevant coherences, (a) 
Non-vanishing off-diagonal elements reaching a steady-state value. Remark that p23 
develops earlier than p4Q. (b) Oscillating coherences that vanish in the steady state 
and therefore do not contribute to the total current. 



two setups: the one we already discussed, where the leads are placed at diagonally 
opposite corners of the sample and a second configuration in which the opposite middle 
sites are used as contacts, like in Fig. 1. The differences are obvious: in the 2nd 
configuration the coupling to the 3rd level decreases and the coupling to the 1st level 
increases. The differences induced in the associated currents are identified in Fig. 9(b): 
The contribution of the 3rd level to the current transmitted in the right lead reduces 
considerably in the 2nd configuration; also the steady-state value of the occupation 
number A^3 is around 0.8 (not shown) which suggests that this level charges more than 
in the 1st setup (A^3 ~ 0.5 in the steady-state). On the other hand the current Jl,i 
increases in the 2nd configuration and the occupation number Ni goes to the steady- 
state value faster that in the previous geometry, meaning that the filling of the 1st level 
is done easier since the opening to the contacts is higher. 
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Figure 9. (Color online) (a) The coupling coefficients |Tg„p in the case of coupling 
at the middle sites and at opposite corner sites. (b) Comparison of the current carried 
by the 1st and 3rd levels from the active region for two configurations of the contacts 
(see the discussions in the text). 
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Figure 10. (Color online) (a) The on-site probability |03(J)P for a 5 x 20 sites sample 
at a magnetic flux $ = 0.2 (b) The same for $ = 0.25. Note the oscillatory behavior 
along the sample, (c) The currents associated to the states shown in (a) and (b). (d) 
The occupation number N3 for the two values of the magnetic field. Other parameters: 
/XL = -1.5 and fiR = -2.25, Vl = Vr = 1, kT = lO"*. 



Another issue we consider in this work is related to the propagation of electrons 
across the sample. More precisely, we want to investigate the behavior of transients 
carried by edge states having similar amplitudes in the contact region but different 
structures along the sample. Such a comparison can be done by changing the magnetic 
flux while keeping the bias window fixed. For example, a 5 x 20 sites plaquette with 
two leads attached to opposite corners has four states located within the bias window 
at $ = 0.2 and $ = 0.25 if the chemical potentials of the leads are fiL = —1.5 and 
Hji = —2.15. The energy of the states within the bias window changes with the magnetic 
field but the tunneling amplitudes {Tgnl"^ remain roughly the same (not shown). The 
main difference occurs in the on-site localization probability |0„(z)p, n = 1,..,4. We 
show in Fig. 10(a) and (b) the localization of the 3rd state from the bias window for the 
two values of the magnetic fiux. At the corners, where the leads are now attached, the 
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states are similar, but the state corresponding to $ = 0.25 has many oscillations along 
the sample edge. The question is then: is there any difference in the transient currents 
carried by the two states? Fig. 10(c) shows that the oscillating state in Fig. 10(b) is 
not favorable for transport, as J/j3 is very small, even if there is a net current entering 
the left lead. In contrast, at $ = 0.2 the 3rd level carries a substantial current in the 
right lead. These results suggest that the state given in Fig. 10(b) should be almost 
filled in the steady-state regime because the electrons are trapped inside the sample 
after entering from the left lead. This fact is confirmed by inspecting the occupation 
number A^3 at the two values of the magnetic flux (see Fig. 10(d)). The steady-state 
value of A^3 ~ 0.5 for $ = 0.2 meaning that on average this level is only half-occupied 
and allows the propagation of electrons in the right lead, whereas at $ = 0.25 the level 
is more filled and therefore the associated current Jr_3 is very small. Note also that for 
$ = 0.25 the steady-state is achieved very slowly. 
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Figure 11. (Color online) The total current Jji for a 5 x 10 sites sample and for 
fast, moderate and slow coupling to the leads. The inset shows the switching functions 
corresponding to the three currents. Other parameters: p(0) = |l)(l|j ^ = 0.2, 
fiL = -1.5 and fin = -2.25, Vl = Vr = 1, fcT = 10-"^. 



We have shown up to now that the transients depend on the position of the 
contacts, on the initial configuration of occupation numbers in the sample and also 
on the localization properties of the carrying state. To complete our analysis we shall 
finally investigate the dependence of the transport properties on the switching functions 
x'. More precisely, we shall keep the form of x' as introduced in the beginning of this 
section and compute the currents for different values of the parameter 7. Fig. 11 gives 
the total current in the right lead for the 5 x 10 sample with p(0) = |1)(1| and with the 
same parameters as in Fig. 4. It is clear that as 7 decreases the filling of the systems 
from the right lead slows down and Jr decreases. For 7 = 0.05, which corresponds to a 
very slow coupling to the leads, we notice a delay even in the filling process (i.e. Jr is 
vanishingly small for t < 20). Also, the small oscillations of the transient that appear at 
7 = 1 are reduced at 7 = 0.25 and disappear as the coupling is established even slower. 
Another interesting fact we could learn from Fig. 11 is that the steady-state current 
does not depend on the switching function. This feature is intuitively expected and was 
even rigorously proved by Cornean et al. [12] 

A natural step forward in our analysis would be to look for memory effects, that is. 
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Figure 12. (Color online) The total currents Jb,l computed from both Markov 
approximation and the non-Markovian versions of GME for a 5 x 25 sites sample. 
Other parameters: p(0) = |1)(1|, $ = 0.2, /i^ = 1.7 and fiR = 1.35, Vl = Vr = 1, 
kT = 10-^ 7 = 0.75. 



to compare the solution of the GME with the Markovian case. Non-Markovian effects 
are expected to appear for strong lead-sample coupling or if one has entanglement in 
the initial configuration. Usually, the memory effects can be neglected if the dynamics 
of the reservoirs is much faster than the dynamics of the sample. In our case we 
cannot check this fact because the relevant time scales (i.e. the correlation time in 
the leads and the dynamics of the sample) are not at hand for the complex systems 
we study. They presumably depend on the bias, on various energy gaps and even on 
the switching functions. This is actually one reason to start from the very beginning 
with a non-Markovian GME. Although a closer investigation of memory effects is left 
for future work we shall give here a brief account on the Markovian version of Eq. (TT3ll 
and a first comparison between the two approaches. Following Timm [10] we replace 
Uo(t,s)p{s)Uoit,sy by p{t) in Eq. ( IT3l) and compute the 11 operator with this ansatz. 
Note that we do not extend the upper integration limit to oo as is done in most papers 
but rather compute numerically the integral. 

We have performed numerical simulations for a 5 x 25 sites system, taking just two 
single particle states within the bias window. The total Markovian and non-Markovian 
currents Jl.r in each lead are shown in Fig. 12. We observe that the memory effects 
appear at t ~ 7.5 and that the currents coincide in the long-time limit t > 45 (the 
switching function reaches its maximum value quite fast, around t ~ 6). Note also that 
the memory effects develop first in J^ and only later in Jj^. This preliminary calculation 
shows that indeed the transient currents may not be suitably captured within the 
Markov approximation, while the steady-state regime is satisfactorily described. Further 
investigation is needed in order to understand memory effects in more complicated 
processes, e.g. pumping. 

Finally we consider as well a two-level system which is extensively studied both 
within GME and rate equation approach. The central region is composed of only 
two sites to which the two leads are attached. The two levels are symmetric w.r.t. 
zero, i.e. Ei = —1,E2 = 1. We have taken a bias window that covers both levels 
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fiL = 1-25, /iij = —1.25. The transients show damped oscillations similar to the ones 
presented by Gurvits et al. [H] In the steady state regime the two states carry equal 
currents. The lowest level oscillates more and its associated current is slightly negative at 
short times due to the back-tunneling processes from the right lead to the sample. This 
feature dissapears when the chemical potential of the right lead is lowered to fiR = —2, 
because the back-tunneling processes are suppresed (not shown). 
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Figure 13. (Color online) The currents in the right for a two-level system. Other 
parameters: p(0) = |1)(1|, $ = 0.0, ml = 1-25 and nn = -1.25, Vl = Vr = 0.25, 
kT = 10-^ 7 = 0.5. 



4. Conclusions and discussion 



One of the main features of mesoscopic systems is that their geometry leaves some 
fingerprints on the transport properties (e.g. the invasive role of the current probes). In 
this paper we have described theoretically the time-dependent transport through such 
structures within the RDO formalism borrowed from quantum optics. When extended 
to open quantum systems this method is a powerful tool for studying electron dynamics 
through and within a sample coupled to biased leads and characterized by a well defined 
initial state. We complement previous approaches [21] |2^ [23l ITT] in the following way: 
i) the GME is solved without the Markov approximation for arbitrary time-dependent 
coupling to the leads; ii) the usual assumption that the spectrum of the sample is entirely 
contained into the bias window is not needed; iii) the transfer Hamiltonian describing 
the coupling between the leads and the sample takes into account the localization of the 
sample states depending on the geometry of the sample and on the region where the 
leads are plugged. 

The GME is solved using the Crank-Nicolson algorithm. In real experiments the 
sample is characterized either by a ground state or by a low-energy excited state and 
for small couplings to the leads one expects that most of the levels below (above) the 
bias window remain occupied (empty) and will not contribute to the current. The 
relevant many-body states are actually few and they are given by all combinations of 
occupation numbers for a bunch of single particle states from the vicinity of the bias 
window. Motivated by this fact one can actually restrict the calculation of the matrix 
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elements for the RDO. A refined version of the argument should hold in the presence of 
the Coulomb interaction as well, using perhaps a Hartree-Fock initial ground state. 

The numerical simulations are obtained for a lattice Hamiltonian. As a main 
application of the method we have computed the transients associated to each level 
of a two-dimensional lattice in the presence of a strong perpendicular magnetic field. 
In this case the currents are carried only by edge states and depend both on the 
contact point and on the topology of the state. Different initial states of the isolated 
system lead to different transients but to the same steady-state current which is not 
achieved at the same time. We have presented a comprehensive analysis of the electron 
dynamics in the transient regime and also studied the relevant matrix elements of 
the RDO (both populations and coherences are discussed). Also we have shown that 
the Markov approximation is appropriate for steady-state calculation but misses some 
memory effects in the transient regime. We want to emphasize that our method not 
only goes beyond the Markov and wide-band approximations, being thus from the very 
beginning more accurate, but it is even more efficient for numerical calculations: the time 
integration can be done recursively, which is not possible in the Markov approximation. 
Consequently the calculations in the Markov approximation took much longer computer 
time than the solution of the GME. 

Some of the transient properties, like the delay in the appearance of a current 
in the drain lead should motivate further experiments in the field. Further expected 
applications of this method include pulse propagation, pumping and time-dependent 
interference effects. 

The electron-electron interaction inside the central region was not included in 
the present calculations and therefore subtle features like charging effects or charging 
effects on the displacement currents could not be discussed here. In the case of many- 
level systems with a specific geometry the treatment of Coulomb interaction within 
a time-dependent framework is a highly non-trivial issue. In particular, the exact 
diagonalization method (mostly used for two-level systems in the cotunneling regime) 
becomes numerically costly due to the large number of many-body states. Therefore 
the calculation of the time- dependent number of particles in the presence of electron- 
electron interactions requires a suitable approximation scheme. We believe however 
that our approach provides a faithful qualitative description of the main dynamical 
processes and gives a first hint about the crucial role played by the spectral properties 
of many-level systems in the transient regime. 

Appendix A. Derivation of GME 

We present here for completeness the derivation of Eq. ( TT3l) following essentially the 
Nakajima-Zwanzig method. To this end we write the equation of motion for W{t) in 
terms of the Liouvillian C: 

^f^dWp_^ ^(t) = ^o + Mi) (A.l) 

dt 
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Co- =[Ho,-], CT{t) = [HT{t),-]. (A.2) 

Next we define two projections: 

P- = PLPRTriTrni-} Q = I - P. (A.3) 

It is straightforward to clieck tlie following properties: 

PCs = CsP, PCT{t)P = 0. (A.4) 

Tlie Liouville equation ( lA.ll) splits then into two equations: 

ihPW{t) = PC{t)PW{t) + PC{t)QW(t) (A.5) 

ihQWit) = QC{t)QW{t) + QC{t)PW{t), (A.6) 

and the second equation can be solved by iterations (T is the time-ordering operator): 

QW{t) = ^ [ dsTexp j--^ f ds'QC{s')\QC{s)PW{s). (A.7) 

Inserting Eq. ( JA.71) in Eq. ( JA.51) and using the properties of P we get the following 
equation: 

ihPW{t) = PCsW{t) 

+^PCT{t)Q I rfsTexpj--^ I ds'QC{s')Q\QCT{s)PW{s). (A.8) 

In order to have an explicit perturbative expansion in powers of the transfer Hamiltonian 
HT{t) one has to factorize the time-ordered exponential as: in the time-ordered 
exponential in Eq. ( JA.81) : 



Texp j-;^ /" ds'QC{s')Q\ = exp{g£oQ}(l + 7^), (A.9) 

where the remainder TZ contains higher powers of Ht- The Born approximation of 
the generalized master equation consists in neglecting TZ. Another technical point is 
that in expanding the unperturbed part exp QCqQ one can replace each Q between two 
Liouvillian Co hj P + Q = 1, since QCqPCqQ = 0. By taking the trace over the leads 
in Eq. ( 1A.8I) one obtains: 

zhp = Csp{t) + ^TiLT^iRiCrit) f dse-'^'-'^'''CT{s)pLPRp{s)}, (A.IO) 

which is nothing but Eq. flT3l) . 
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